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Abstract 

The qf t++ package is a library of C++ classes that facilitate numerical (not algebraic) quantum field theory calcu- 
lations. Mathematical objects such as matrices, tensors, Dirac spinors, polarization and orbital angular momentum 
tensors, etc. are represented as C++ objects in qf t++. The package permits construction of code which closely re- 
sembles quantum field theory expressions, allowing for quick and reliable calculations. 



1. Introduction 

It is often desirable to describe particle physics 
processes using a covariant tensor formalism. This 
formalism contains Dirac spinors and matrices, po- 
larization and orbital angular momentum tensors, 
etc. For cases involving higher spin particles or large 
values of orbital angular momentum, the presence 
of high rank tensors can make explicit calculation of 
the desired quantities (typically scattering or decay 
amplitudes) cumbersome at best and impossible at 
worst. 

A common approach is to perform algebraic 
manipulations to reduce the number of tensor 
contractions prior to coding the expressions of in- 
terest. There are several software packages avail- 
able to facilitate this approach, e.g. FeynCalc [1] 
and FeynArts [2]. Packages like GRACE [3] and 
CompHEP [4] take this a step further by (mostly au- 
tomatically) producing distributions from physical 
models. While these packages are extremely use- 
ful in certain areas of physics {e.g. writing event 
generators), they arc not ideal for performing an 
event-based partial wave analysis (PWA). 

In this type of analysis, numeric values for ampli- 
tudes must be calculated for upwards of 100 million 
events (data and/or Monte Carlo). Thus, the am- 
plitude calculation software must not only be easy 
to use, but also be capable of performing numerical 



calculations of expressions involving matrices and 

high-rank tensors using as little cpu time as possible. 
The qf t++ package satisfies both of these criteria. 

The operations of matrix multiplication and ten- 
sor contraction (the building blocks of all such calcu- 
lations) can be broken down into nested loops; thus, 
they are easy to perform (numerically) on a com- 
puter. The qf t++ package was developed to take ad- 
vantage of this fact by performing numerical calcu- 
lations (not algebraic manipulations), making these 
types of calculations easier and more accessible to a 
larger portion of the physics community. The qf t++ 
package has been used in a number of PWA's to date 
(see, e.g., [5]) and has performed exceptionally well. 

One further point before discussing specific details 
about the software, this package was designed to 
calculate tree- level expressions. No work has been 
done to provide a simple way of performing loop 
integrals. 

2. Overview 

The the main design goal for the application pro- 
gramming interface (API) was to make the code re- 
semble the mathematical expressions as closely as 
possible. To facilitate this goal, each type of mathe- 
matical object {e.g. tensors, spinors, etc.) has a cor- 
responding C++ object in qf t++. Through the use 
of operator overloading, operations such as matrix 
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multiplication and tensor contraction are handled 
by the objects themselves; not by the user, i.e. the 
user never has to keep track of any indices. 

As a very simple example, consider the electro- 
magnetic vertex eu{pi,'mi)"f'^u(p2,m2)ef^(p~^,m~^). 
After variable declaration and initialization (dis- 
cussed below), this expression can be written using 
qf t++ as follows: 



e*Bar (ul (ml) ) *gaiiima*u2 (m2) *eps (mg) ; 



The object types determine how "multiplication", 
i.e. the * operator, is to be performed. The matrix 
multiplications and tensor contractions are handled 
internally by the objects. This feature allows for 
self-documenting code, greatly reducing the proba- 
bility for mistakes. 



3. Basic Operations 

The two main types of operations required in 
quantum field theory calculations, but not con- 
tained in standard C-|— 1-, involve matrices and 
tensors. From these constituents, it is easy to build 
classes which handle Dirac matrices, covariant pro- 
jection operators, etc. This section describes how 
the qf t++ package implements these basic types of 
operations. 

This package makes heavy use of template classes; 
thus, a number of template utilities have been devel- 
oped for performance and/or API reasons. Among 
these are compile-time detection of inheritance and 
parameter passing optimization [6], along with se- 
lective inclusion of class methods [7] . 

Unlike many standard C-l— I- template classes, all 
of the template classes defined in the qf t++ pack- 
age are designed such that instantiations of differ- 
ent types are fully compatible, provided the types 
themselves have the nec;essary operators defined, 
e.g. the following code is legal: 



SomeClass<Tl> scl; 
SomeClass<T2> sc2; 
scl*sc2; 



if the operator * is defined between types Tl and T2. 



3.1. Matrix Operations 

Matrix operations are handled by the template 
class Matrix, which can store any data type which 
can be stored in a CH — h STL container class, i.e. 
it can store any object which can be stored by 
std: : vector. A complete set of methods are pro- 
vided including those useful in Quantum Field the- 
ory calculations, such as Trace and Adjoint, along 
with all the necessary operators. 



3.2. Tensor Operations 

Tensor operations are handled by the template 
class Tensor, which can also store any data type 
which can be stored in a C-|— I- STL container class; 
however, the most useful types are typically double 
and complex<double>. Methods are provided to 
perform Lorentz transformations, symmetrization 
and a number of other functions. 

A number of operators are provided to perform 
tensor contractions. The Tensor class assumes that 
all indices are either raised or lowered. The field 
theory objects discussed in Section 4 use the for- 
mer, i.e. they are contravariant. The * operator per- 
forms standard multiplication if either object is not 
a Tensor and contraction of a single index other- 
wise. For example, if x is a rank-2 Tensor, then the 
code snippet 3*x simply multiplies each of the ele- 
ments of x by three; however, the snippet x*y will 
evaluate the expression x^^^y^P if y is also a rank- 
2 Tensor. Full contraction of all possible indices is 
performed using the I operator; thus, x | y evaluates 
the expression x^^y^^. 



Symbolic Expression 


qft++ Code 




x*y 




x*/.y 


X ^ 2//ll/l2---^n 


(xly) 



Table 1 

Example tensor operations using the Tensor template class. 



Additional classes are provided for the Minkowski 
metric, (MetricTensor), and the Levi-Civita 
tensor, ep,uaf3 (LeviCivitaTensor). The Vector4 
template class, which is derived from Tensor, pro- 
vides a number of additional methods specific to 4- 
vectors, e.g. CosTheta and Beta. 
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4. Object-Oriented Field Theory 

From the Tensor and Matrix classes discussed 
above, all of the necessary objects for quantum field 
theory calculations can be constructed. In this sec- 
tion, a brief overview of the formalism will be pre- 
sented followed by a discussion of the corresponding 
qft++ class. More detailed descriptions of the for- 
malism are given in [8,9]. Recent examples of apply- 
ing this formalism to PWA - for which this code was 
developed - can be found in [5,10,11]. 

4.1. Integral Spin Wave Functions 

The wave function of a particle with integral spin- 
J, 4-momentum p and spin projection to some quan- 
tization axis m, is described by a rank-J tensor, 
^tii...ii.j{p,'>n)- The Rarita-Schwinger conditions for 
integral spin- J are 



p^'€f,,f,^.„f,^...i^j{p,m) = 



■^IJ.l...lM--f' 



(la) 
.nj{p,m) (lb) 



(p, m) = 0, 



(Ic) 



for any /ij , jij , and reduce the number of independent 
elements from 4"^ to (2 J -|- 1). 

The spin- J projection operator, defined as 



piJ) 



jn 

is used to construct the particle's propagator. 

As a simple example, consider a massive spin-1 
particle. The Rarita-Schwinger conditions for spin-1 
simply require p'' c ^,(p,m) = 0. Thus, in the parti- 
cle's rest frame the energy component of the wave 
function is zero. The spatial components are then 
chosen to be 

ei;±l)=T^(l,±t,0), ei;0) = (0,0,1). (3) 

The wave function can be obtained in any other 
frame through the use of Lorentz transformations. 
The spin-1 projection operator is given by 



spin is set using the constructor, e.g PolVector 
eps(3) would be used for a spin-3 particle. The 
PolVector object must then be initialized for a 
given 4-momentum. At this point, the user can de- 
cide whether or not the particle is to be on-shell. For 
example, if a spin-1 particle is initialized on-shell 
then the projection operator is given by 



Pil\p) 



PtJ.P,y 

p2 ' 



(5) 



otherwise, it is calculated using (4). This option is 
also available in qft++ for half-integral spin parti- 
cles. 

After initialization, the sub-states can be ac- 
cessed via calls like eps (m) , which returns the 
Tensor<complex<double> > object for sub-state 
TO. The spin-J projection operator can be ac- 
cessed easily in the code using the method 
eps.ProjectorO, which returns an object of type 
Tensor<coinplex<double> > with a rank of 2 J . 

4.2. Half-Integral Spin Wave Functions 

The wave fimction for a spin- 1/2 particle is de- 
scribed by a 4-component Dirac spinor, denoted 
as u{p,m), where p and m again represent the 4- 
momentum and spin projection of the particle. The 
representation chosen here leads to the following 
form of the spinors: 



u{p, m) = \/E + 




(6) 



where E{w) is the energy(mass) of the particle and 
x(to) are the standard non-relativistic 2-component 
spinors. 

Wave functions for paxticles with higher half- 
integral spin, denoted as u^^„,^j_^^^{p,m), are con- 
structed using tensor products of integral spin wave 
functions and the spin- 1/2 spinors described above. 
The Rarita-Schwinger conditions for half-integral 
spin-J are [8] 



{'y'^Pf,-w)u^^. 



(p, to) = 



= u, 



IJ,l...lMj...lJ,i...lXj. 



-1/2 fe™) 



where w is the mass of the particle. 

Wave functions for particles with integral spin 
are handled in qf t++ by the PolVector class. The 



p'''u^,,...^,...^j_^^^{p,m)=Q 
'y'''u^^...^,,...u.J_,/AP^rn) =0 



(7a) 

(7b) 
(7c) 
(7d) 
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qft++ Clahb 


Symliol 


( 'oilCC'pl 


Matrix<T> 


dij 


matrices of any dimension 


Tensor<T> 




tensors of any rank 


MetricTensor 




Minkowski metric 


LeviCivitaTensor 




totally anti-symmetric Levi-Civita tensor 


DiracSpinor 


"a'i-.-;^j-i/2(P' 


half-integral spin wave functions 


DiracAntiSpinor 


v{p, m) 


spin- 1/2 anti-particle wave functions 


DlracGamma 






DiracGammaS 


7^ 


Dirac matrices 


DiracSigma 






PolVector 


eMl...M,7 (P.H 


integral spin wave functions 


OrbitalTensor 




orbital angular momentum tensors 



Table 2 

A partial list of qft++ 



, a complete list can be found at [12]. 

4.3. Dirac Matrices 



ff''*^^UMl...M....K,.-M,7-l/2(P>^) =0, (7c) 

for any /z j , Hj , and reduce the number of independent 
elements from 4-^+^2 to (2 J + 1). 

The spin- J projection operator is then defined as 



P^-^^ (p) 



1 

2w 



X! "''I i^j-in{P^ rn)u^^ {p, m). (8) 



Classes are also provided to handle the Dirac 
matrices, 7'' (DiracGainma) and 7^ (DiracGeomnaS), 
along with a^"" = |[7'',7''] (DiracSigma). Each 
of these classes is derived from the common base 
class Matrix<Tensor<complex<double> > >. Thus, 
they inherit all of the necessary Matrix and Tensor 
operators. 



For example, the spin-1/2 projection operator is 

simply P(^/2)(p) — 2^(p''7^+ui), while for spin-3/2 
the projection operator is given by 

X (p^IHp) + lP^'J{p)rPl,l\ph''') . (9) 

Wave functions for particles with half-integral 
spin are handled in qft++ by the DiracSpinor 
class. The spin is set using the constructor, 
e.g. DiracSpinor u(3/2.) would be used for 
a spin-3/2 particle. As with the PolVector 
class, the DiracSpinor class must be initialized 
for a given 4^momentum. The sub-states can 
then be accessed via u(m), which returns the 
Matrix<Tensor<complex<double> > > object for 
sub-state m. The spin- J projection operators can be 
easily accessed using the method u.ProjectorO 
which returns a 4 x 4 Matrix of rank-(2J — 1) 
Tensor<complex<double> > objects. We also note 
here that the quantity u{p, m) is obtained using the 
function Bar (u (m) ) . 



4.4. Orbital Angular Momentum Tensors 

Two particles, with 4- momenta pa and p^. can be 
coupled to a state of pure orbital angular momen- 
tum, using the operators i|L^i^2-- w The total and 
relative momenta are defined as P = pa + Ph and 
Pab = \ {Pa — Pb) respectively. The orbital-angular- 
momentum operators are then built using the rela- 
tive momentum and the spin-^ projection operator 
as follows: 

^\illi2.--lie ^liiH2..-neviV2...ve{P)PabPab ■ ' -Pab' i^^) 

These operators satisfy the Rarita-Schwinger condi- 
tions 







■^IJ,iH2...m...Hj--.ne /:ii([i2--.(Uj •••w 



1X2-. -Hi. .-Hi... III. 



0, 



(11a) 
(lib) 
(11c) 



for any /Xj, jij, which insure that they have {2^ + 1) 
independent elements. 
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In the qf t++ package, orbital-angular-momentum 
operators are handled by the OrbitalTensor class. 
This class inherits from Tensor<double>; thus, 
after construction it can be used just like any 
other tensor. Setting the tensor elements of -L/j/...^^ 
(OrbitalTensor object orbL) for 4-momenta Pa 
and ph (Vector4<double> objects pa and pb) is 
done by simply calling orbL.SetP4(pa,pb). 

4.5. Additional Utilities 

Functions are also provided to calculate use- 
ful quantities such as Clebsch-Gordon coefficients, 
Wigner D-functions, Breit-Wigner and Regge prop- 
agators, etc. 

5. Exeimple Applications 

In this section, a few simple examples will be ex- 
amined. The qf t++ package computes the values of 
expressions numerically; thus, the 4-momenta of the 
particles involved must be known. In an event-based 
partial wave analysis, these would be obtained from 
the experimental data and/or Monte Carlo events. 
In the case that one wants to calculate theoreti- 
cal angular distributions, cross sections, polarization 
obsc^rvables, etc.. events with the desired kinemat- 
ics must be generated as input for the qf t++ code. 
In the examples below, the assumption is made that 
one of these methods is employed. For the example 
plots shown in this section, the latter method was 
used. 

5.1. X{2-) ^ujK^ tt+w-tt°K 

Consider the decay of a particle, X, with spin- 
parity = 2~ into an oj and K via F-wa.ve. The 
invariant decay amplitude for this process is propor- 
tional to 

A oc e*(p,,m^)i(3)^-«(p^^)e,„(P,M), (12) 

where Puj,muj{P,M) are the momentum and spin 
projection of the uj{X) and p^^K is the relative mo- 
mentum of the (jjK system. 

For this example, assume that the X is produced 
via e+e~ annihilation resulting in population of only 
the M = ±1 sub-states. The decay distribution is 
then obtained by calculating the intensity 

lex ^ J2 (13) 

M=±lm„=±l,0 



To calculate this distribution using qft++, the 
necessary variables must first be declared: 



PolVector epso; // omega 
PolVector epsx(2) ; // X 
OrbitalTensor orb3(3) ; // L~3 
Tensor<coinplex<double> > amp; 
Vector4<double> p4o,p4k,p4x; 



For each point at which J is to be calculated, the 
4^momenta must be set and the polarization states 
initialized using calls to SetP4; however, if the kine- 
matics are such that the X mass is constant, then 
p4x and epsx will only need to be initialized once. 

In the code, a loop would then be performed over 
all values of cos 9 (the decay angle of the w in the 
X rest frame) for which the decay intensity is to 
be calculated (or over events). At each point, the 
calculation would be performed as: 



double intensity = . ; 
for (Spin m = -1 ; m <= 1 ; m+=2){ 
f or(Spin mo = -1 ; mo <= 1 ; mo++){ 
amp = conj (epso(mo))*orb3|epsx(m) ; 
intensity += norm(ainp() ) ; 

> 

} 



In this way, the value of I can be calculated at any 
number of points in cos^ (or, for any number of 
events). 

It is worth examining this more carefully. The 
tensor contractions in the code above are performed 
by first evaluating the * operator which contracts 
^*^{Pui,'m'w) into the first index of L^^^^'"^{puik)- 
The result is a rank-2 Tensor whose two in- 
dices are contracted via the I operator into 
both indices of e^a{P,M). The result is a rank-0 
Tensor<complex<double> > whose value is ac- 
cessed via the () operator. 

To check that the code is working for this simple 
example, the intensity in (13) can be calculated in 
the X rest frame by making a slight modification 
to the non-relativistic helicity formalism solution as 
follows: 

I(x^|/(i;^/u;„,A)(3OlA|2A)cii,;,(0)|2, (14) 

M,X 



5 



I — analytic 




i. r 

-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 

cos(e) 



Fig. 1. (Color Online) Intensity vs cos 8: Angular distri- 
bution calculated for the decay X{2~) — > wK using qft++ 
(black filled squares) compared with an analytic solution 
(red line). See text for details. 

where Ei^,Wuj,X are the energy, mass and hehcity of 
the uj and /(x, ±1) = 1, f{x,0) = x accounts for 
the effects of the boosts on the covariant w hehcity 
states. Notice that as w^, the non-rclativistic 

solution is recovered. 

The expression in (14) can then be rewritten 
purely in terms of the decay angle as follows: 

I(x (2 cos^ 61- 1)2 +cos^ 61 

+9^— ^ sin^flcos^^. (15) 
J 

Figure 1 shows the angular distributions obtained 
by calculating (12) using qft++ compared to the 
modified helicity formalism expression given in (15) , 
normalized to have the same integral as the qf t++ 
solution. The two calculations give the same angular 
distributions, i.e. the code is working properly. 

This example can be extended by considering the 
secondary decay uj —^ Tr+Tr^Tr*^. The amplitude for 
this process can be written as 

-4^^7r+7r-7rO « ie^^„^p^+p"_ p^o e'" (Pc^ , m^) , (16) 

which, in the lu rest frame, simplifies to 

This is the standard non-relativistic result [13]. 
To incorporate this decay, (12) must be rewritten 

as 

A « ujML^^^'"""iP^K)eUP,M), (18) 
where. 



« „2 „,„2 , r ^^pcppUp.-p^o, (19) 

and Wi^^T^ are the mass and width of the u). 

To add this decay to the qft++ calculation, the 
following variables need to be defined: 



complex<double> i(0,l) ; 
Teiisor<complex<double> > omega; 
LeviCivitaTensor levi; 
Vector4<double> p4pip , p4pim , p4pi0 ; 



Then, for each point the intensity is calculated as: 



omega = epso.Projector()*levi 

*p4pip*p4pim*p4piO 

*BreitWigner (p4o , . 78256 , . 00844) ; 
double intensity = 0. ; 
for (Spin m = -1 ; m <= 1 ; m+=2) { 

amp = omega*orb3 I epsx(m) ; 

intensity += normCampO) ; 

} 



Since all of the objects are covariant, no extra 
boosts or rotations to the lo rest frame are required. 

5.2. 7rp A — > 7rp 

As an example involving half-integral spin par- 
ticles, consider the reaction up — > A(1232) np. 
The invariant scattering amplitude for this process 
is proportional to 

M cx uipf,mf)p''fPl3\p)py{p„m,), (20) 

where Pi{pf) and mi(mf) are the 4- momentum and 
canonical spin projection of the initial(final) proton 
respectively. For this example, the mass-dependence 
of the propagator will be ignored. In the code, adding 
this would simply involve multiplying the amplitude 
by a complex<double>. 

The distribution in the scattering angle, i.e. 
the angle between the initial and final protons 
(cos 9 = Pi ■ pf), is obtained by calculating the scat- 
tering intensity 

I(x ^ \M\^ (21) 

mi,mf 
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To calculate this distribution using qft++, the 
following variables must be declared: 



DiracSpinor ui ,uf ; // proton spinors 
DiracSpinor delta(3/2 . ) ; 
Matrix<complex<double> > £mip(l,l) ; 
Vector4<double> p4_i ,p4_f ,p4_delta; 



For each point at which X is to be calculated, the 
4-momenta must be set and the spinors initialized 
using calls to SetP4; however, if the 4-momentum of 
any given particle does change from point to point, 
then its corresponding object would only need to 
be initialized once. 

In the code, a loop would then be performed over 
all values of cos d for which the scattering intensity 
is to be calculated (or, again, over all events). At 
each point, the calculation would be performed as: 



double intensity = 0. ; 

for(Spinm_i = -l/2. ; m_i<=l/2. ; m_i++){ 
for(Spinm_f = -1/2. ; m_f <= 1/2. ; m_f++)-[ 
amp = Bar(uf (m_f ))*p4_f 

♦delta . Pro j ector () *p4_i*ui (m_i) ; 
intensity += norm(cmip(0,0)) ; 

} 

} 



In this way, the value of X can be calculated at any 
number of points in cos B. 

If the A (1232) projector is on-shell (see discussion 
in Section 4.1), then the numerical calculations can 
be checked in the overall center-of-mass frame using 
the non-relativistic helicity formalism. In this frame, 
the scattering intensity is proportional to 



MIa,(^)I 



^ oc 1 + 3 cos^ I 



(22) 



where Ai(A/) are the initial (final) proton helicities. 

Figure 2 shows the angular distributions obtained 
by calculating (20) in the overall center-of-mass 
frame using qf t++ compared to the non-relativistic 
helicity formalism calculation used to obtain (22). 
The result obtained using (22) was normalized to 
have the same integral as the covariant calculation. 
Clearly the two calculations give the same angular 
distributions. 
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■ analytic 
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Fig. 2. (Color Online) Intensity vs cos 6: Angular distribu- 
tion calculated for the reaction 7rp — > A(1232) — > 7rp using 
qf t++ (black filled squares) compared with an analytic solu- 
tion (red line). See text for details. 



5.3. Compton Scattering 

As a final example, the unpolarized Compton 
scattering cross section will be calculated to leading 
order in a. The two well-known tree-level ampli- 
tudes, corresponding to s- and u-channel electron 
exchange diagrams, are 



A. 



"f''e^{k,mj)u{p,mp) (23a) 



i) — t' + w 



where p{p') and k(k') denote the initial(final) elec- 
tron and photon momenta respectively and w is the 
mass of the electron. The full scattering amplitude, 
from which the cross section can be calculated, is 
obtained by combining these two processes. 

The calculation of this scattering amplitude is 
similar to that of the previous example. First, the 
following variables must be declared: 
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DiracSpinor ui ,uf ; // e spinors 
DiracSpinor u_ex; // exchanged e spinor 
PolVector epsi , epsf ; // photon pol .vecs 
DiracGamma gamma; // gamma'mu 
// s- and u-chcLnnel propagators 
Matrix<complex<double> > prop_s; 
Matrix<complex<double> > prop_u; 
Vector4<double> pi ,pf ,ki ,kf ; // 4-momenta 



A loop would then be performed over all values of 
cos 9 (scattering angle in the lab frame) for which 
the cross section is to be calculated. At each point, 
the propagators are obtained using the following 
code: 



u_ex.SetP4(pi+ki, 0.511) ; 
prop_s = u_ex. Propagator ; 
u_ex.SetP4(pi-kf ,0.511) ; 
prop_u = u_ex. Propagator ; 



The scattering intensity is then obtained by looping 
over spin projections, Spin m_ef ,m_ei ,m_gf ,m_gi, 
and calculating the amplitudes given in (23) as 



ainp_s = Bar(uf (m_ef ))*gamma 

*conj (epsf (m_gf ) ) *prop_s*gamma 

*epsi (m_gi) *ui (m_ei) ; 
ainp_u = Bar(uf (m_ef ))*gamma*epsi(m_gi) 

*prop_u*gamma*conj (epsf (m_gf ) ) 

*ui (m_ei) ; 



To get the cross sections the appropriate scale fac- 
tors (a^, phase space factors, etc.) must then be 
applied. 

Figure 3 shows the differential cross section for 
a 10 MeV incident photon calculated using qf t++ 
compared to the well-known spin-averaged Klein- 
Nishina formula [14]. Both methods give the same 
results. 

6. Discussion 

The examples in this paper were chosen because 
they can also easily be solved analytically, allow- 
ing for comparisons of the two calculations. In gen- 
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cos(e) 



Fig. 3. (Color Online) ^^"^ ^ (barns) vs cos 9: Compton scat- 
tering differential cross section for a 10 MeV incident photon 
calculated using qf t++ (black filled squares) compared with 
the Klein-Nishina equation (red line). See text for details. 

eral, this is not the case; however, even for very 
complicated expressions involving high rank tensors 
or symmetrization, the qf t++ calculations are still 
very manageable (more complicated examples can 
be found online [12]). The code is also optimized for 
performance. The examples given in this paper run 
at ~ 3 kHz, i.e. on a 2.5 GHz processor, one can 
compute results at about 3000 kinematic points per 
second. This package makes it possible for any physi- 
cist to compare and/or fit a theoretical model to his 
or her data by coding up the expressions themselves. 

It is clear from these examples how this package 
could be useful for a partial wave analysis, but it 
can also be used for other types of analyses. For ex- 
ample, it can be used to perform numeric checks of 
analytic calculations. The qft++ package can also 
be used for cases where analytic solutions are desir- 
able but not feasible. For example, consider the case 
where an expansion of an expression in powers of 
some variable, x, is needed but an analytic solution 
is not available. The qf t++ package could be used to 
evaluate the expression at a large number of values 
of X. The coefficients in the expansion could then be 
extracted by fitting the qf t++ results. 

The qft++ package is available to all users from 
http : / /www-meg . phys . emu . edu/qf t| f+. 
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